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We characterize a tree's spatial foliage distribution by the local leaf area density. 
Considering this spatially continuous variable allows to describe the spatiotemporal 
evolution of the tree crown by means of 3D partial differential equations. These offer a 
framework to rigorously take locally and adaptively acting effects into account, notably 
the growth toward light. Biomass production through photosynthesis and the allocation 
to foliage and wood are readily included in this model framework. The system of 
equations stands out due to its inherent dynamic property of self-organization and 
spontaneous adaptation, generating complex behavior from even only a few parameters. 
The density-based approach yields spatially structured tree crowns without relying on 
detailed geometry. We present the methodological fundamentals of such a modeling 
approach and discuss further prospects and applications. 

Keywords: leaf area density, continuity equation, functional-structural plant model, crown plasticity, competition 
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1. INTRODUCTION 

In terms of model scale, light sensitive functional-structural 
tree growth modeling has experienced the emergence of various 
trends. Organ-level approaches bring about a high precision of 
physiological processes, averting inaccuracies and effects of scale 
non-invariance, which may arise from simplifications in larger- 
scale approaches. Moreover, arbitrary small-scale biophysical or 
biochemical processes can in principle be readily induced. The 
LIGNUM model (Perttunen et al., 1996; Sievanen et al., 2008), 
the model by Sterck et al. (2005) or the L-Peach model (Allen 
et al., 2005) are examples of this model category. Local light inter- 
ception (cf. also Chelle and Andrieu, 2007 for a methods review) 
determines the production and allocation of biomass. Their detail 
of physiological and morphological processes is at the same time 
the drawback of these models. On the one hand, the large num- 
ber of organs implies high computational costs — all the more in 
competition scenarios of multiple trees. On the other hand, their 
detail can make these models susceptible to the propagation of 
errors, which could have been compensated for in an averaging 
rough scale approach. 

Other organ-level models like Greenlab (Yan et al, 2004; 
Cournede et al., 2006) a priori focus on the topology in terms 
of the plant's structure. This implies the inability of easily tak- 
ing physiological and structural responses to varying local light 
conditions into account. As one consequence, the approach can- 
not be straightforwardly applied to scenarios of competition for 
light: An additional yet non local competition index provides for 
this (Cournede et al, 2008). In another context of formal gram- 
mars used for tree growth simulation, Kurth and Sloboda (1999) 
present the 2D concept of the shadow-relevant cone of a shoot in 
order to take local light conditions into account, which in turn 
affect the rewriting system. Comparable methods have been used 
by Purves et al. (2007) and Takenaka (1994). 



Sonntag (1996) presents a model in which the spatial motion 
and allocation of leaf area is based on heuristic rules on a 2D 
cellular space. Though quite different in terms of formalism, his 
approach bears conceptual resemblances to ours. 

Models with a rougher scale use to impose certain charac- 
teristics on the crown shape. For instance, the Balance model 
Grote and Pretzsch (2002) and the model by Sorrensen-Cothern 
et al. (1993) describe the crown shape in terms of disk-like hor- 
izontal layers. This technique implies advantages with regard to 
the computational speed as well as a general robustness, com- 
pared to small-scale models. Yet it does so at the expense of a 
thorough plastic spatial crown structure. When applied to compe- 
tition scenarios, these models often make use of empirically fitted 
competition indices (e.g., Pretzsch, 1992). 

While being affiliated more with the latter forestry mod- 
els in the attempt to describe crown structure and dynamics 
macroscopically for applications at the stand level, the present 
approach attempts a middle course between the fine organ- 
centered way and the rougher pre-imposition of the crown 
form. We characterize a tree's spatial foliage distribution via 
its local leaf area density. The focusing on this variable cir- 
cumvents difficulties in terms of robustness in the geometri- 
cally detailed models accounting for individual leaf positions. 
At the same time, the locality allows for arbitrary spatial 
structures. 

Applying Beer Lambert's law allows to express the local light 
conditions within a crown as a function of the local leaf area 
density. Aiming at an increased future light interception, local 
leaf area density is assumed to tend to move toward the light. 
This approach, which induces the spatial expansion of the crown, 
translates directly into a partial differential equation. Details are 
specified based on mass conservation and optimization consid- 
erations. The technique notably allows to account for a local and 
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spontaneous adaptiveness with regard to changing environmental 
light conditions. 

The mathematical approach of density-based partial differen- 
tial equations has long established itself in spatial biology and 
ecology (cf. e.g., Okubo and Levin, 2002). In the context of 
macroscopic individual plant modeling, so far notably in the form 
of diffusion equations, it has proven applicable for root growth 
and proliferation (see Page and Gerwitz, 1974 for the original 
approach, Reddy and Pachepsky, 2001 for a review of later devel- 
opments and Dupuy et al., 2010 for a current advance). A 2D 
diffusion approach for the foliage of crops with the objective to 
model competition in different field densities, without consider- 
ing the vertical dimension, was presented by Beyer et al. (2014). 
Partial differential equations can generate, even from only a few 
terms, complex self-adapting dynamics, which is indeed present 
in biological systems. Attempts to reproduce this by simpler terms 
often requires a larger model framework and set of parameters. 
The latter, in turn, requires large data sets, which are often not 
available. 

In this article we will exemplify the use of the leaf area 
density-based partial differential equation approach by means of 
a simplified model setup presented in section 2. The continuity 
of the approach suggests embedding into in the context of con- 
tinuously growing trees sensu (Halle et al, 1978), i.e., "with no 
marked endogenous cessation of extension [and] a more or less 
constant production of leaves and/or shoots throughout the year" 
(Barthelemy and Caraglio, 2007). 

Since merely selected dynamics are presented, without con- 
structing a realistic model, we settle for illustrations of key qual- 
itative properties of the approach instead of quantitative data 
comparison. Throughout the article, we will point out possibili- 
ties of extending and customizing process assumptions while pre- 
serving the advantages of the overall methodological framework. 
Extending future steps are discussed in section 3. 

2. MODEL FRAMEWORK 

For clarity's sake all recurrently appearing variables and param- 
eters and their definitions are listed in Table 1. We describe the 
spatial foliage distribution of a tree canopy by means of the 



Table 1 | Key model variables and parameters. 



ff(x, t) 


Leaf area density in x at time f (in m 2 m~ 3 ) 


a(; t) 


The map K 3 -* R, x i-> a(x, f) 


A 


Mean leaf transmittance 


Mx, v) 


Cosine of the angle between leaf plane normal and sun ray 


si 


Upper half unit sphere jx s R 3 :|jx|| = 1 , X3 > 0) 




Energetic efficiency (in gMJ -1 ) 


PAR (v, t) 


Radiation from direction v at time f (in MJ m~ 2 s~ 1 ) 


b(x, t) 


Local biomass production in x at time t (in g s~ 1 ) 


SLA 


Specific leaf area (in m 2 g~ 1 ) 


WD 


Wood density (in g m~ 3 ) 


P 


Pipe model theory constant (in m 2 m~ 2 ): 1 unit a = P Units 




pipe cross-sectional area 


Mr 


Length of the sapwood pipe leading to x (in m) 


L(x, f) 


Local radiation in x induced by the leaf density a(t) (in 




MJrrr 2 ) 



(local) leaf area density a(x, t) > 0 (in m 2 m~ 3 ) in a point x = 
(x\ , X2, X3) € R 3 (each entry in m) at a time t > 0 (in years), i.e., 
the spatial density of the total one-sided green leaf area (previ- 
ously considered e.g., by Sinoquet et al., 2005). The map a(-, f) : 
R 3 — > R, x h-> a{x, t) is continuous for any t. For brevity, we will 
use the notion leaf density in place of leaf area density. 

We aim to describe the evolution t \-+ a{-, t). To this end we 
will first determine the biomass production B of a tree, which, 
along with the senescence S of old biomass, allows to describe the 
net biomass increment of a tree corresponding to a leaf density 
a(-, t) at a given time t: 

|-m(t) = B(t) - S(t) 
at 

This, in turn, will be distributed among foliage and sapwood 
according to the pipe model theory (Shinozaki et al, 1964), 
specified in the subsequent paragraph. The coaction of foliage 
allocation and senescence in the case of a continuously growing 
tree induces what we abstractly interpret as a continuous motion 
of a(-, t), in particular directed toward the light, aiming at an 
increased future biomass production. This perspective leads to the 
description of the course of 1 1— > a(-, t) essentially by means of a 
continuity equation. 

Sapwood Associated to a Leaf Density: The pipe model theory by 
Shinozaki et al. (1964) allows to determine the sapwood mass cor- 
responding to an arbitrary leaf density a (•, f). In the present con- 
text, the theory states that for any point x € R 3 with a (x, t) > 0, 
a sapwood pipe, in charge of the transport of water and nutrients, 
leads from x down to the roots with a length denoted by -r , its 
cross-sectional area being proportional to the leaf density a (x, t) 
at its tip via a constant P. The mass of the pipe leading to x then 
equals 

a(x, t) ■ P ■ \\x\\y ■ wd, 

cross-sectional area length 

WD denoting the wood density (in gm~ 3 ). 

It is worth mentioning that, motivated by the limitations 
to the pipe model theory, pointed out e.g., by Tyree (1988), 
Pouderoux et al. (2001), and Deleuze and Houllier (2002), gen- 
eralizations have been suggested: A noteworthy approach is the 
one by Bouchon et al. (1997), who, based on an allocation per- 
spective, reason that the pipe does not necessarily have a constant 
cross-sectional area along its path, but more generally a one 
exponentially decreasing toward the stem base. This principle 
integrates in our context with only minor technical changes. 
Likewise, the approach by Letort et al. (2008), parametrically 
combining the pipe model approach and a uniform, common 
pool sapwood allocation, would be feasible. 

As for the pipe's length ||x|| x from x to the root tip, for the sake 
of simplicity, here we follow the multi-species approach used by 
Sonntag (1996), based on the branch architecture of coniferous 
species and the assumption that root length equals branch length, 
resulting in 

\\x\\ r :=\x 3 \+2-\\(x l ,x 2 )\\ . (1) 
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A more accurate choice for a particular species can be made by 
taking specific characteristics of its branching geometry (such as 
branching angles) and topology into account. 

Finally, the sum of foliage and sapwood mass is given by 



f a(x, t) f 

I — dx+ I a(x, t) ■ P ■ \\x\\r ■ WDdx 

Jk 3 SLA J R 3 



(2) 



folias 



apwood i 



where SLA denotes the specific leaf area (in m 2 g *). 
2.1. BIOMASS PRODUCTION 

We determine the amount of biomass produced through 
photosynthetic activity by a given leaf density. To this, we 
take direct and diffuse radiation into account, cf. Fu and 
Rich (1999), using the horizontal celestial coordinate system, 
with R 2 x {0} being the local horizon of a tree rooting in 
(0, 0, 0), and the vector (1, 0, 0) pointing north. The unit 



directional vector er(f) € Si_ := {x G 



1, x 3 > 0}, 



under which the sun is seen from the tree at daytime t reads 
a(t) = (cos( - Az(f)) • cos(Alt(f)), sin ( - Az(f)(f)) • cos(Alt 
(f)), sin ( Alt (t))), where Alt(f) and Az(f) denote the time 
dependent altitude and azimuth, respectively. 

For diffusive radiation, a uniform diffuse model (uniform 
overcast sky) is applied, in which incoming diffuse radiation is 
assumed to be the same from all sky directions. Let PAR^;,- (f) 
and PARjiff (f) (in MJm" 2 ) denote the photosynthetically active 
direct and diffuse radiation at time f, respectively. Then the total 
radiation from direction v € S 2 , at time f is 



PAR (v, t) := PARdiff (f) + PAR dlr (f) • 1 



{v=cr(f)| 



(f), 



with the indicator function 1^ (x) := 1 if x e A and 0 else. 
2.1.1. Isolated Tree 

Incoming radiation is partly intercepted by the tree's foliage and 
partly passes through it. The fraction between 0 and 100% of radi- 
ation from direction veSj which actually reaches the point x e 
M 3 can be determined using Beer-Lambert's law, where foliage 
characterized by leaf density acts as a light absorbing medium 
with locally varying a -concentration. This fraction reads 



exp I — A 



/ 



(3) 



where the extinction coefficient A < 1 represents the mean light 
transmittance of foliage (Monteith, 1969; Nouvellon et al., 2000) 
and 

X(x, v) := N(x) ■ v 

takes into account the angle between the sun ray and foliage in x, 
N(x) G S 2 ^ denoting the unit normal to the plane in which foliage 
in x lies. N(x) can be chosen according to leaf angle distribution 
models without further ado (Wang et al., 2007) provide a review. 

Assuming that local biomass production is proportional to the 
total locally intercepted ratiation via an energy efficiency [i (in 



gMJ -1 ) (Monteith, 1977), the instantaneous local biomass pro- 
duction b(x, f) (in gm~ 3 s~ : ) in x € R 3 at time t for the leaf 
density a (•, f) reads 



b(x, t) = fi ■ 



v) ■ a(x, t) ■ exp 



(- A 7 



MS, v) ■ a(S , f) d§ ■ PAR (v, t) dv. (4) 



radiation reaching .v from direction v 



radiation from direction v intercepted in x 

2. 1.2. Population of Trees 

At a time f, let ai(-, ?),..., a n (-, t) denote the leaf densities of 
n trees, shading each other and competing for light. Then the 
instantaneous local biomass production of tree i e {1, . . . , n] in 
x generalizes to 



bi(x, t) = [i 



/I 



(Xi(x, v) • <Xi(x, t)) 2 



exp 



The fraction 



\ x+R + -v-' =1 



•PAR(v, t)dv 



™ — —. — r — -. — r is the part of the incoming radia- 

2_,j=i kj{x,v)-otj{x,t) 

tion in x that is attributed to tree i's foliage in x. In particular it 
reduces to 1 if two trees' crowns do not occupy common space. 



2.2. DYNAMICS 
2.2.1. Mass balance 

We determine the instantaneous change in living mass (conduc- 
tive sapwood and foliage mass) due to the production of new, 
and the senescence of old biomass. For convenience we assume 
that the senescence of leaves as well as the loss of conductivity of 
sapwood depend on time only. If sapwood and foliage that have 
existed for Tw and rp years become nonconductive and senescent, 
respectively, then the living mass at time t reduces by 



S(t) 



f 



- t ce(x, t - rp) 



SLA 

P ■ \\x\\r ■ WDdx. 



dx + 



a(x, t — Tw) ■ 



At the same time it increases by the total instantaneous biomass 
production at f, i.e., B(t ) := L 3 b(x, t) dx. Thus we have 



dt 



m(t) = 5(f) - S(t) 



(5) 



for the change in living mass at time f. A priori, the possibility 
i(t) < 0 is not excluded for arbitrary parameters. However, 



St 1 



since m(t) > S{t), it follows that m(t) > 0 and thus B(t) > 0 for 
all t > 0. 
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2.2.2. Leaf density dynamics 

With this information on the global mass of the tree at hand, we 
consider the variable 



a(x, t) 



a(x, t) 



,3 a(x, t)dx 



i.e., the leaf density modulo mass, or (mass-)relative leaf density, 
with the property that / R3 a (x, t ) dx = 1 for all f , which is under- 
stood as an indicator of the spatial structure of the real leaf density 
a(-, t). Instead of describing the course of a(-, f) directly, we do so 
for (?(•, t) and deduce a(-, t) as a(-, t) = k(t) ■ a(-, t), where X(t) 
is chosen such that the living mass corresponding to k(t) ■ «(•, t) 
(cf. (2)) equals indeed m(t) yielded by (5). Thus 



Ut) 



m(t) 



■Ik 



dx + / R3 a{x, t) ■ P ■ 



\X\\X ■ WD 



dx 



(6) 



The continuity of the growth process of the trees we consider 
suggests a description of the course of a(-, t) in terms of a 
(mass-conserving) continuity equation, 



3* 



a(x, t) = Vjc • (p(x, t), 



(7) 



in which the relative leaf density is subject to a transport motion 
induced by a continuous flux <p : R 3 x R + — > R 3 , itself deter- 
mined by a{-, t). The idea of this transporting flux <p is that it 
incorporates both the effect of the allocation of new leaves and of 
the abscission of old ones on the spatial structure of foliage, which 
in combination, induces what we describe as an abstract motion 
of leaf density. 

A predominant driver in the spatial dispersal of leaf density 
is the local expansion toward the light, aiming at an increase in 
future light interception. We formally embed this factor in the 
above framework, where it will take on the role of the flux </>. For 
some given leaf density a(-, t) let L : R 3 x R + — > R be defined 
by 

L(x, t)= I X{x, v) ■ exp ( - A • / X(? , v) • ot(|, t) d$\ ■ 

J Si \ Jx+R + -v ) 



The term V x • <f> describing a movement of leaf density toward 
the light contains spatial derivates of a function of integrals over 
a, which makes (7) a partial integro- differential equation. 

2.3. SIMULATIONS 

In this section we illustrate some structural properties of the 
model. For convenience we simplified radiation to be vertical 
only. Some details of the numerical implementation of the model 
are presented in the appendix. 

Figures 1, 2 illustrate the evolution of leaf density in the course 
of time, as well as the vector field <f>. The term is sensitive to any 
change in the local light conditions induced by shading; (p instan- 
taneously adapts and points in the direction of the greatest light 
increase. Aside the general spatial expansion toward the light, we 
notably observe the predominant presence of foliage at the crown 
hull rather than its interior. 

2.3.1. Population of Trees 

Together with the adjustments in terms of biomass production 
addressed at the end of section 2.1, the approach generalizes to 
competition scenarios when a(-, f) is replaced by a i('> f ) m 
(8), alongside the different initial states ai(-, 0), ot2(-, 0), . . .. We 
illustrate the dynamic effects to which competition gives rise by 
means of a simplified scenario: Consider a sufficiently large stand, 
in which the trees' stem bases, as a point set in R 2 x {0}, generate 
a Voronoi-tesselation which is regular. If radiation is assumed to 
be radially symmetric (as done e.g., by Perttunen et al., 1998), 
the analysis of all competing trees reduces to that of a single 
one for which periodic boundary conditions for (7) are added 
on the boundary of the tree's 3D cell, i.e., the extension of the 
appropriate 2D Voronoi cell in the x$ -dimension. 




I r r s 



PAR (v, t) dv 



(8) 



0 



FIGURE 1 | (x q . x 3 )-cross section of a leaf density-characterized crown 
in the course of time, subject to the present model dynamics, 
corresponding to t = 10 (A), 20 (B), 30 (C) and 40 (D) years with 
simulation parameters adopted from Letort et al. (2008). A darker color 
indicates a higher leaf density a. The arrows in (D) indicate the vector field tj>. 



The function L measures the intercepted light in x per m 2 leaf 
area. The gradient V X L points locally in the direction of the great- 
est rate of increase of intercepted light. In addition, similar to 
Beyer et al. (2014) we define the flux to correspond to the existing 
leaf density in x, so that finally we have 



(j>{x, f) = k ■ a(x, t) ■ V x L(x, t) 



(9) 



for a mobility constant k. This local gradient approach is moti- 
vated by "the observation that a tree is capable of acquiring 
also gradient information about its environment and that growth 
might be directed along these gradients (Schmidt and Wulff, 
1993; Aphalo and Ballare, 1995)" (Sonntag, 1996). 























FIGURE 2 | 3D view of Figure 1D. 
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Periodic boundary conditions induce that the light conditions 
on the other side of the boundary are considered identical to those 
within, accounting for another tree growing in equal measure and 
shading its environment. 

This implies that, when, for simplicity, further assuming 
the essential light incidence to be vertical, periodic boundary 
conditions reduce to no flux conditions on the boundary: There, 
due to the identical light conditions on the other side of the 
boundary, the light gradient V X L, governing the flux <p, changes 
from pointing further outwards to a zero flux. 

Figure 3 shows the different stages of this scenario for an 
underlying square tesselation. 

The overly sharp boundary between two crowns is a conse- 
quence of the unrestrained mobility of foliage in this simplified 
approach, and would fade when additional features correspond- 
ing to a branch structure are included, cf. section 3. 

The conspicuous concentration of leaf density at the upper 
edges of the canopy in the competition case (cf. Figure 3C) results 
from the two factors of (i) the regular expansion of leaf density 
uninfluenced by the boundary condition and (ii) the immigration 
of leaf density from lower regions whose horizontal expansion 
had abruptly turned into a vertical one after reaching the Voronoi 
cell boundary. We are unable to provide evidence or counterev- 
idence to determine whether this phenomenon corresponds to 
reality. In any case, as observable in Figure 3D, this is only an 
intermediate state, before in the long-term a homogeneously dis- 
tributed high concentration of leaf density in the top canopy 
establishes itself. 

Wheras an isolated tree grows radially symmetric, this sym- 
metry is eventually broken for a non-circular Voronoi cell, such 
as the square used just now. Figure 4 illustrates this. 

3. DISCUSSION AND FUTURE PROSPECTS 

The aim of this article was to show how local leaf area den- 
sity, a concept opposing geometrically detailed individual leaf 
configurations, can be used to approach macroscopic tree crown 



dynamics. Its integration into the formal framework of partial 
differential equations allowed to rigorously formulate the growth 
toward light. In the simulations we observed the generation of 
self-organization and adaptiveness that come along with this 
modeling approach. The simplistic model framework was meant 
to draw attention to the key mechanisms and their dynamic 
effects. 

Foliage dynamics are by nature coupled to the tree's branch 
structure, which has not been taken into account on a topo- 
logical or geometric level in this article. Future work, with the 
aim of introducing more spatial heterogeneity to the approach, 
begins here. Taking merely the stem and the most vigorous pri- 
mary branches into account while leaving the finer structures to 
the leaf area density concept may already suffice to tackle crown 
plasticity satisfactorily. While in the present simplified model the 
motion of local leaf area density is governed by light only and oth- 
erwise unrestrained, a simple branch architecture can add sort 
of directional inertia to that motion, channeling local leaf area 
density in several major directions, resulting in heterogeneous 
foliage clustering, representing individual branches. In particu- 
lar, this includes the incorporation of genetically predetermined 
branching angle spectra. Introducing a branch structure, even 
if only a rough one, is moreover accompanied by a refinement 
of the pipe length term (1), governing the distribution of mass 




FIGURE 4 | Cross-section of (A) a diagonal plane (B) the x 3 )-plane 
at f = 15 years. (C) 3D view. 







FIGURE 3 | (xi, x 3 )-cross section and 3D view of the stages in a competition scenario at t = 10 (A), 20 (B), 30 (C) and 40 (D) years with periodic 
bounday conditions at 3.5 m. Due to properties of self-symmetry, visible in Figure 1, these stages are in fact qualitatively similar no matter the cell size. 
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between foliage and wood according to the pipe model theory, 
thus determining secondary growth. 

Taking into account the organization of growth units, i.e., as 
weakening our assumptions on neoformation and polycyclism in 
section 2.1.1, or considering immediate vs. delayed bud outbreak, 
would bring the model closer to actual tree architecture dynamics. 

Alongside the phototropism considered in this article, more 
biomechanical constraints that have feedback influences on the 
growth, such as hydraulic aspects sensu (Ryan and Yoder, 1997; 
Tyree and Zimmermann, 2002), in particular in the context of 
growth limitation, the avoidance of interlocked growth due to 
mechanical stress of touching branches (Oliver and Larson, 1990) 
or gravitropism represent perspectives for model extensions. 

The application and validation of a refined model based on the 
theoretical framework presented in this paper will benefit from 
empirical data on local leaf area density. Conceivable ways to 
obtain this include the following three, which are currently being 
practically explored: Firstly, from the direct recordings of local 
light intensities at various positions [x- 1 ', . . . , x^ n >] C K 3 within 
a canopy, the map 

{x m , ...,x (n) } R> 0 
x i — > a(x) 

for a discrete, but arbitrary fine domain can be obtained by 
applying Beer-Lambert's law in the reverse way. Secondly, high- 
definition multi-directional 3D terrestrial laser scan data and 
appropriate skeletonization algorithms allow to relate to a leafless 
tree a set of cylinders representing branch segments (Raumonen 
et al., 2013). Automatically removing all scan points correspond- 
ing to cylinder diameters (i.e., branch thicknesses) above a certain 
threshold allows to spatially isolate recent shoots, for which then a 
relation to foliage can be assumed. Thirdly, scanning a tree at the 
beginning of spring during the process of bud opening, when — 
in contrast to the case of a fully developed canopy — laser rays still 
reasonably penetrate the crown, and deducting from this image 
the point cloud yielded by a scan of the completely leafless tree, 
may allow to obtain a local bud density, from which the the leaf 
density can be deduced. Thus, assessing the spatial foliage distri- 
bution in terms of local leaf area density in a functional-structural 
crown dynamics model can make good use of this data type for 
model calibration and validation. 

The property of locally and spontaneously adapting to chang- 
ing light conditions suggests that the present partial differ- 
ential equation approach can be applied to competition sce- 
narios both in pure as in mixed tree groups. Empirical find- 
ings, based on laser scans, about the plasticity Bayer et al. 
(2013) may thus be approached from a functional-structural 
modeling point of view. These perspectives are currently being 
explored. 
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4. APPENDIX 

NUMERICAL IMPLEMENTATION 

A finite volume scheme Toro (2009) is used, in which we con- 
sider Oyi- := /j.j a ( x > t)dx (with a similar notation for the 
other variables) on a regular mesh with cells E,* = [iAx, (i + 
l)Ax[x[jAx, (/' + l)Ax[x[kAx, (k+ l)Ax] for i,j,keZ and 
Ax « 1. 

The local light incidence (8) (and similarly the local biomass 
production (4)) for vertical radiation, v = e$, is computed as 

Lijk(t) = Xijkiei) ■ exp I — A ■ Ax ■ ^Xj/fcte) • Uijk(t) ■ 



PAR(e 3 , t)dv 



The general case for v € is conceptually similar, yet technically 
more extensive: For a given v, the above sum reaches over all cells 
EjW whose intersection with the line through the center of the 
cell EjW and pointing in the direction v is non-empty, taking 



into accout the individual length li>jw(v) < Ax of this intersec- 
tion by means of the additional coefficient in the sum 
term. 

As for the PDE, let 



01,yt( f ) = k 



oii+\jk{t) + otijk(t) L i+l j k (t) - Lijk(t) 



Ax 



and 4>2,ijk' $3,»Tfc be defined accordingly for the other spatial 
entries. The standard discretization of the divergence is given by 

div^(f) = i : i + 



+ 



Ax 
Ax 



Ax 



and finally, using the Euler method for Af <SC 1, we obtain 

Oi/jfc(f + Af) = a y7 t(f) + Af • divyjt(f). 
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